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Abstract 

A genetic algorithm approach suitable for solving multi-objective optimization problems is described and 
evaluated using a series of simple model problems. Several new features including a binning selection 
algorithm and a gene-space transformation procedure are included. The genetic algorithm is suitable for 
finding pareto optimal solutions in search spaces that are defined by any number of genes and that 
contain any number of local extrema. Results indicate that the genetic algorithm optimization approach is 
flexible in application and extremely reliable, providing optimal results for all optimization problems 
attempted. The binning algorithm generally provides pareto front quality enhancements and moderate 
convergence efficiency improvements for most of the model problems. The gene-space transformation 
procedure provides a large convergence efficiency enhancement for problems with non-convoluted 
pareto fronts and a degradation in efficiency for problems with convoluted pareto fronts. The most difficult 
problems — multi-mode search spaces with a large number of genes and convoluted pareto fronts — 
require a large number of function evaluations for GA convergence, but always converge. 


Nomenclature 


F 

F* 

f 

f 

G n 

IR 

ISELECT 

ITRAN 

NSEED 


N c 

Nco 

N g 

Nm 

N 0 

P 

Pi 

P2 

R 

R(0,1) 


set of scalar objective functions 
optimal set of F values 
scalar objective function 
maximum value of f 
n th GA generation 
integer ranking array 

user-specified integer parameter controlling which selection algorithm is used 
user-specified integer parameter controlling the gene-space transformation option 
user-specified parameter that controls how many solutions with different initializing 
random number generator seeds are averaged together to construct a single 
convergence history curve 

fixed number of chromosomes in each GA generation 

number of constraints 

number of genes in each chromosome 

number of modes (hills or peaks) associated with the MP1 and MP2 model problems 
number of scalar objective functions 

user-specified vector with four elements that controls which modification operators are 
used in going from G n to G n+1 

user-specified parameter controlling the probability that a specific gene will be modified 
using the perturbation mutation operator (0 <p-, <1) 

user-specified parameter controlling the probability that a specific gene will be modified 
using the global mutation operator (0 < p 2 < 1) 
basis vector matrix with elements 

random number generator that returns a random value between 0 and 1 
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u 


■^max, 

P 


unitary transformation matrix with elements Ujj 
i th gene from the chromosome from the n th GA generation 
j fl chromosome from the n th GA generation 
user-specified maximum limit on the i th gene 
user-specified minimum limit on the i th gene 

user-specified parameter controlling the size of the perturbation mutations (0 < J3< 1) 


subscripts 

/ 

j 

k 

m 


gene or decision variable index 

chromosome index 

objective function index 

mode index (associated with MP1 and MP2) 


superscripts 

n GA generation or population index 

t temporary chromosome and gene values obtained after selection but before operator 

modification 


Background 

Numerical methods for optimizing the performance of engineering problems have been studied for many 
years. Perhaps the most widely used general approach involves the computation of sensitivity gradients. 
These methods — called gradient methods — have been utilized to produce optimal engineering 
performance in a wide variety of different forms. The reliability and success of gradient methods generally 
requires a smooth design space and the existence of only a single global extremum or an initial guess 
close enough to the global extremum that will ensure proper convergence. 

In contrast to gradient-based methods, design space search methods such as genetic algorithms (GA) — 
also referred to as evolutionary algorithms (EA) — offer an alternative approach with several attractive 
features. The basic idea associated with the GA approach is to search for optimal solutions using an 
analogy to the theory of evolution. The problem to be optimized is parameterized into a set of decision 
variables or genes. Each set of genes that fully defines one design is called an individual or a 
chromosome. A set of chromosomes is called a population or a generation. Each complete design or 
chromosome is evaluated using a “biological-like” fitness function that determines survivability of that 
particular chromosome. For example, in aerospace applications, the genes may be a series of geometric 
parameters associated with an aerospace vehicle that is to be optimized for payload delivered to orbit, 
aerodynamic performance or structural weight. The fitness function takes as input all the geometric 
parameters and returns a numerical value for the fitness — the size of the payload, the aerodynamic 
performance or the structural weight. 

During solution advance (or “evolution” using GA terminology) each chromosome is ranked according to 
its fitness vector — one fitness value for each objective. The higher-ranking chromosomes are selected to 
continue to the next generation while the lower-ranking chromosomes are not selected at all. The newly 
selected chromosomes in the next generation are manipulated using various operators (combination, 
crossover or mutation) to create the final set of chromosomes for the new generation. These 
chromosomes are then evaluated for fitness and the process continues — iterating from generation to 
generation — until a suitable level of convergence is obtained or until a specified number of generations 
has been completed. 

Constraints can easily be included in the GA optimization approach either by direct inclusion into the 
fitness function definition — a so-called penalty constraint — or, more commonly, by including one or more 
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constraints into the fitness function evaluation procedure. For example, if a design violates a constraint, 
its fitness is set to zero, i.e., it does not survive to the next generation. Because GA optimization is not a 
gradient-based optimization technique, it does not need sensitivity derivatives. It theoretically works well 
in non-smooth design spaces containing several or perhaps many local extrema. 

General GA details including descriptions of basic genetic algorithm concepts can be found in Goldberg , 1 
Davis , 2 and Beasley, et al. Additional useful studies, which survey recent activities in the area of genetic 
algorithm or evolutionary algorithm research including the presentation of model problems useful for 
evaluating GA performance are given in Deb , 5 Van Veldhuizen and l_amont 6 and Jimenez, et al . 7 

A disadvantage of the GA approach is expense. In general, the number of function evaluations required 
for a GA optimization process to converge, exceeds the number of function evaluations required by a 
finite-difference-based gradient optimization (see the results presented in Obayashi and Tsukahara 8 and 
Bock 9 ). This situation is offset, to an extent, by the ease with which GAs can be implemented in parallel or 
distributed computing environments. 

Despite being relatively new, genetic algorithms have already been applied in many applications over a 
broad range of engineering fields. A brief survey of single discipline applications is presented in Holst and 
Pulliam , 10 and will not be discussed further in this report. 

Other applications involving GA search methods have been made in the area of multi-objective or multi- 
discipline optimization, i.e., optimization problems in which two or more objectives are simultaneously 
optimized. These methods, referred to as MOGA (multi-objective genetic algorithm) methods, are 
especially attractive because they offer the ability to directly compute so-called “pareto optimal sets” in a 
single computation instead of the limited single design point traditionally provided by other methods. 

The pareto optimal set or pareto front, as it is common called, includes optimal solutions for each of the 
individual objectives, as well as a range of tradeoff solutions in between, which are themselves optimal 
solutions. Providing a range of solutions to a multi-objective optimization problem is a powerful approach 
because it allows the designer to see the effect of decision variable variation on the design space in the 
form of optimal tradeoffs. Thus, the designer can choose individual objective weighting factors after their 
full influence is quantitatively known. 

In recognition of the importance of the MOGA approach, many theoretical developments have been 
published in recent years. In particular, Van Veldhuizen and Lamont , 6 and Zitzler, et al . 11 present 
reasonably complete surveys of MOGA methodology with a special emphasis on how to compare GA 
performance from one algorithm variation to another. Other researchers, including Deb , 5 Deb, et al ., 12 
Van Veldhuizen 13 and Lohn, et al ., 14 have developed and/or utilize a suite of test problems as a standard 
in evaluating MOGA convergence efficiency as well as the accuracy of the final pareto front. A key aspect 
of pareto front development is diversity. Does a particular MOGA produce good coverage over the entire 
pareto front or are some regions poorly resolved while other regions have high levels of undesirable 
clustering? This issue is studied in Laumanns, et al . 15 and Deb, et al . 16 

Still other MOGA issues are associated with archive strategies. As the pareto front develops, many 
solutions are found that lie on the pareto front that cannot be retained within the fixed population size of 
many schemes. How to retain or archive this information for the benefit of the MOGA while maintaining 
reasonable bounds on the archive file has been studied in Knowles and Crone . 17-18 One last example of 
MOGA research lies in the area of an interactive algorithm development. Parmee, et al . 19 present a 
MOGA approach that allows changes to be made in GA operators as wells as in objective definitions as 
the MOGA advances from generation to generation. 

One area of MOGA research that is even more voluminous than the area of theoretical developments is 
that associated with GA multi-objective applications. Examples of multi-objective optimization applications 
include airfoil optimization by Marco, et al ., 20 Naujoks, et al ., 21 Quagliarella and Della Cioppa , 22 Vicini and 
Quagliarella , 23 Hamalainen, et al . 24 and Epstein and Peigin , 25 missile aerodynamic shape optimization by 
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Anderson, et al., 26 wing optimization by Anderson and Gebert, 27 Sasaki, et al., 28 Oyama, 29 Obayashi, et 
al., 30 and Ng, et al. 31 

Additional MOGA examples in the area of turbomachinery optimization include rocket engine turbopump 
design by Oyama and Liou 32 and compressor blade design by Benni 33 and Oyama and Liou. 34 In some of 
these examples the multiple objectives were obtained by considering two different aerodynamic design 
points. In others the multiple objectives involved different disciplines including aerodynamics, structures, 
controls and/or electromagnetics. 

One last area of MOGA research that bears mention is in the area of hybrid methods, the utilization of 
MOGA optimization in conjunction with another type of optimizer. This is an attractive area of research 
because MOGA methods are particularly good at finding a global extrema, but not particularly good in 
converging the optimal solution to tight levels of accuracy. Briefly, several examples where hybrid 
approaches have been utilized include Giotis, et al. 35 and Tursi 36 where a MOGA has been coupled to a 
neural network for aerodynamic shape optimization, and Vicini and Quagliarella 37 and Brown and Smith 38 
where a MOGA has been coupled with a gradient optimizer. 

Various definitions and the multi-objective genetic algorithm used in the present study are described next. 
Details associated with each of the operators, including selection, passthrough, random average 
crossover, perturbation mutation and mutation are presented. MOGA convergence efficiency is then 
evaluated using several model problems from two general points of view — the effect of design space 
characteristics on GA convergence and the effect of GA control parameter specification on GA 
convergence. 


Problem Statement: Single Objective Optimization 

A single-objective optimization problem can be stated as follows: Let f be a scalar function of N G 
independent variables, x, , defined on some domain Q 

f=f(X)=f(x u -,x h -,x Ng ) (la) 

In this notation X is the vector of design space decision variables. The maximum value of f , indicated by 
f * , is obtained by finding the values of X = X* such that 4 

f* = max {f } = f(X*) = f(x ; , • • • ,x* , • • • ,x„ G ) (1b) 

The above maximization operation is subject to N co conditions or constraints indicated by 

c„(X)<0 n = \2,---,N co (1c) 

The constraints placed on the decision variable vector X by Eqs. (1c) essentially serve to limit the design 
space within Q for which the optimal solution is sought. 


Problem Statement: Multi-Objective Optimization 

For optimization problems involving more than one objective, which are simultaneously optimized, the 
situation is more difficult. This is because each objective must play a role in determining the optimal 
solution. In the optimization process, conflicts might arise among the various objective functions, i.e., the 


4 For the purpose of simplifying the discussion of algorithmic details, maximization is generally assumed. 
The logic for minimization is a straightforward modification and will not be discussed. 
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optimal values of each individual objective, in general, will not occur for the same decision variable 
vector, X . As a result, the “optimal solution” for a multi-objective optimization problem is typically a range 
or a set of solutions, which represent tradeoffs in objective space. 

To determine when one solution is better than another for multi-objective problems the concept of 
dominance is utilized. 1 A vector U = 11(0,, ■■■,u N ) is said to dominate another vector 
V = V(i/V • • y, v • • y w ) if and only if u, > v, for all / and there exists at least one value of / such that u ,■ > v, . 
A vector defined on some domain Q that is not dominated by any other vector defined on Q is said to be 

non-dominated on Q . 

A multi-objective optimization problem can be stated as follows: Let F be a set of N 0 scalar objective 
functions, f k , each dependent upon the same decision variable vector, X , which is defined on some 
design space Q 


F=F (fj(X),...,f k (X),-,f Ab (X)) (2a) 

As above, the decision variable vector X consists of N G independent components. The multi-objective 
optimization problem involves finding the set of X =X* that produce non-dominated values for F =F* on 
Q . This set of values F* is called the pareto optimal set or the pareto front. 

For each F the constraints 


c n (X)<0 n=\2,-,N co (2b) 

must be satisfied. Existence of these constraints serves to limit or reduce the size of Q for which the 
optimal solution is sought. 


Genetic Algorithm 

The genetic algorithm optimization procedure utilized to solve the multi-objective optimization problem, as 
described by Eqs. (2), is now presented. It is closed related to the GA optimization procedure presented 
in Holst and Pulliam, 10 which was designed for single objective problems. As mentioned in the 
introduction, the general idea behind GA optimization is to discretely describe the design space using a 
number of decision variables, x, . In GA parlance these parameters are called genes, and the / subscript 
refers to the gene number. Each set of genes that leads to the complete specification of an individual 
design, i.e., each decision variable vector, X , is called a chromosome and is indicated by 


Xj = Xy(X 1 ",-",X,”v 


'<Xn 


j) 


( 3 ) 


where the j subscript, added to X , identifies the chromosome number. In addition, the j subscript has 
been added to each gene, so as to indicate which chromosome each gene value is identified with. The n 
superscript had been added to indicate the GA generation number, which is iteratively advanced as the 
solution converges. Thus, X" is the y th chromosome for the n th generation that consists of N G genes. 


For aerodynamic shape optimization problems, the design space genes are typically a series of 
geometric parameters, e.g., airfoil thickness and camber and/or wing sweep, twist and taper. For many 
GA applications genes are computationally represented using bit strings and the operators used to 
manipulate them are designed to accommodate bit string data. In the present approach, following the 
arguments of Oyama, 39 Houck, et al. 40 and Michalewicz, 41 real-number encoding is used to represent all 
genes. The key reason for using real number encoding is that it has been shown to be more efficient in 
terms of CPU time relative to binary encoded GAs. 41 In addition, real numbers are used for all genes in 
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the present implementation because many engineering applications involve decision variables that are 
best described using real numbers, e.g., the geometric parameters in aerodynamic shape optimization. 
Thus, using real number encoding eliminates the need for binary-to-real number conversions. 

Initialization 


Once the design space has been defined in terms of a set of real-number genes, the next step is to form 
an initial generation, G°, represented by 


G° = (X 1 °,-,X°,..,X° c ) (4) 

where N c is the total number of chromosomes. Each gene within each chromosome is assigned an initial 
numerical value using a process that randomly chooses numbers between fixed user-specified limits. For 
example, the / th gene in an arbitrary chromosome is initialized using 

Xj = R{0, 1)(X maXf — -^min, ) + ^min, (5) 

where x max . and x min are the upper and lower limits for the /' th gene, respectively, and R(0,1) is a 
random number generator that delivers an arbitrary numerical value between 0.0 and 1.0. 

The random number generator used in the present study requires an integer input — a seed value. If the 
integer is positive, the next number in the current random number sequence is always returned. If the 
integer is negative, the random number sequence is reset. Utilization of the same negative seed value 
will always reset the random number generator to the same random sequence. Each new solution begins 
by resetting the random number generator using a single call to R(0,1) with a negative seed value. All 
other calls to R(0,1) during that solution use a positive seed value. Thus, a solution can be repeated by 
simply using the same initial seed value or rerun to determine statistical variation by using a different 
initial seed value. 

Fitness evaluation 


After a generation is established — either the initial generation or any of the succeeding generations in the 
evolutionary process — fitness values, F ; n , are computed for each chromosome using a suitable function 

evaluation. This is analogous to the objective function evaluation in gradient methods and is represented 
using 


F? = F(X") ( 6 ) 

For example, for a multi-discipline optimization (MDO) problems involving the simultaneous maximization 
of two separate and distinct objectives, ^ and f 2 , the fitness evaluation represented by Eq. (6) consists of 
the following 


fr=wj) 

f2=f 2 (Xj) 

where, for example, the first function ^ might be the aerodynamic drag of an aerospace vehicle 
(constrained to fly at fixed lift) and the second function f 2 might be the structural mass of the same 
vehicle. Once all the genes in a specific chromosome have been specified, f\ can be evaluated using a 
suitable CFD solver to obtain the drag and f 2 can be evaluated using a suitable structural analysis routine 


6 



to obtain the structural mass. In this case, of course, the optimization problem would be one of 
minimization. 

Ranking 

The purpose of the ranking operation is to determine a set of integer values for the ranking array, IR . 
One integer value is required for each chromosome. The ranking algorithm is quite different for multi- 
objective optimization problems relative to its single-objective counterpart. As such, both situations will be 
discussed in this section. The ranking array values — once determined — are then used in the GA selection 
process. 

Single Objective Optimization — The GA ranking algorithm for single objective optimization problems is 
quite simple. Whichever chromosome has the highest fitness is ranked number one (IR = 1), whichever 
has the second highest fitness is ranked number two (IR = 2), and so on. This ranking algorithm, more 
formally stated, is given by 


ic = 1 

if(f j n <f l n )ic=ic + 1 / = 1, -,A/ C 
IRj 1 = ic 


j = 1, • Wc 


( 7 ) 


where j and / are special counters that range over all N c chromosomes in the current population or 
generation level. 

Multi-Objective Optimization — For multi-objective optimization problems the ranking procedure is more 
complex and utilizes the concept of dominance, which was defined previously. 

A chromosome with a fitness vector F that is not dominated by any other fitness vector within the design 
space is said to be a non-dominated chromosome. The optimal solution set F* or pareto front includes all 
solutions with fitness vectors that are non-dominated and that satisfy the constraints given by Eqs. (2b). 
The numerical approximation to the pareto front must involve a suitably complete set of discrete solutions 
so as to describe the optimal values of each individual objective — these are the pareto front endpoints — 
as well as the non-dominated tradeoff solutions in between the individual objective’s optimal values. 

The use of dominance for determining pareto fronts in multi-objective optimization problems can be 
clarified by considering a simple example involving two-objectives without constraints. When N a = 2 the 
optimization problem given by Eq. (2a) can be restated as 

F = F(/j(X)f 2 (X)) (8) 

The solution set for the problem defined in Eq. (8) contains all solutions with fitness vectors F = F(/j/ 2 ) in 
the feasible domain that are non-dominated — or for numerical optimization, a suitably populated set of 
discrete solutions with fitness vectors that are non-dominated. For a design point (^V 2 *) to be non- 
dominated, there can exist no points in the design space (/j/ 2 ) such that 

/j > /j * and f 2 >f 2 * 
or 

and f 2 >f 2 * 

This situation is presented schematically in Fig. 1 where the shaded region is the feasible region of the 
design space and the thick dark line contains all members of the pareto front. The square symbols 
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represent the maximum values of each of the individual objectives and also serve to define the pareto 
front endpoints. 



*• 


INFEASIBLE 

REGION 


PARFTD PROMT 
max fi 


Fig. 1 Schematic of a two-objective design space in which ^ and f 2 are simultaneously maximized. The 
thick dark line is the pareto front. 

With the concept of dominance in hand the actual ranking process can now be presented. There are a 
number of ranking procedures available for use in multi-objective optimization. The one used in the 
present study is called Goldberg ranking. 1 After all of the fitness values have been determined, each 
chromosome is checked for dominance. Those chromosomes that are non-dominated are given a 
number one ranking (IR = 1) and then set aside. The remaining chromosomes that are non-dominated 
are given a number two ranking (IR = 2) and then set aside. This process continues until each 
chromosome has an integer value for the ranking array, IR . In general, with this approach, the number of 
different integer values contained within the ranking array will be small, at least small in comparison to 
N c , as many chromosomes will be ranked near the top with a value of 1 , 2 or 3. 

The above procedure, by itself, represents a legitimate algorithm for the ranking operation. However, 
there is a refinement that provides a considerable enhancement to GA convergence. Before describing 
this refinement, two additional concepts must be defined — the active file and the accumulation file. 

The active file is simply a specific name used for the current generation of chromosomes, that is 

G n =(X?,-,x;,...,X^) 

is the n th generation active file for the GA iteration process. As the GA solution evolves, the active file is 
always fixed in size at N c chromosomes. The first N 0 elements in the active file — for the present 
implementation — always contain the chromosomes that posses the maximum fitness values for each of 
the N 0 objectives. 

The accumulation file is the list of all non-dominated chromosomes that have been found from all 
generations combined. As the GA solution evolves, the accumulation file typically grows in size with more 
chromosomes being added after each new generation. As new non-dominated chromosomes are added 
to the accumulation file, old chromosomes that lose their dominance must be deleted, thus ensuring that 
the accumulation file contains nothing but number-one ranked chromosomes. Because each 
chromosome stored in the n th generation accumulation file is non-dominated, it serves to define the 
pareto front or, at least, the n lh generation approximation to the pareto front. The use of an accumulation 
file makes sense only when N 0 >2. 
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The multi-objective ranking procedure described above utilizes only the chromosomes in the active file, 
i.e., each chromosome in the active file is ranked relative to the other chromosomes in the active file. But 
this philosophy potentially wastes a plethora of information because not every number-one ranked 
chromosome for multi-objective optimization problems can be retained from one generation to the next in 
the active file. 

This is where the refinement in the ranking routine enters in. Once each chromosome in the active file is 
ranked using the standard routine, an additional test is performed to see if any number-one ranked 
chromosomes in the active file are dominated by any of the chromosomes in the accumulation file. If this 
is the case, then the ranking number associated with the newly ranked chromosome in the active file is 
decremented by one. This ensures that the ranking routine produces number one rankings that are 
number one in a global sense. 

Selection 


After the ranking array is established, with or without the accumulation file option, the GA algorithm 
passes to the selection process to determine which chromosomes will continue on to the (n + 1) st 
generation and which will not. In the present study three different selection operators will be studied: (1) a 
new and simple technique called greedy selection, (2) a traditional algorithm called tournament selection, 
and (3) a third algorithm based on a “binning” strategy called bin selection. In each case the selection 
process picks established chromosomes — either from the n th generation active file, X", or from the 
accumulation file — and then presents them to the modification process, which will be described shortly. 

Greedy selection — This selection algorithm is quite simple. It selects all of its chromosomes from the n lh 
generation active file, i.e., from X" . It is implemented using the following 


/= 1 


if (IRj < / ) then 
X]=X n j 
1=1 + 1 


\J = \N C 


\it = \N c 


endif J 

if (/ > N c ) stop 


( 9 ) 


where each selected chromosome X] is placed in a temporary holding array indicated by 

G f = (X' v .-,X',-X c ) 

Note how the highest ranked individuals in the r? th generation are selected multiple times, individuals with 
average ranking are selected a small number of times, and individuals with the lowest rankings are not 
selected at all. This biasing toward individuals with the highest rankings is a key element in any GA. The 
chromosomes represented by G f are next used by succeeding modification operators to produce G n+1 . 

Tournament selection — The tournament selection algorithm is quite simple and is used widely — one 
variation or another — in many GA optimization applications. The present variation is implemented as 
follows. Using a random process, three chromosomes are chosen from the n th generation active file, i.e., 
from X” . The ranking array values, IR , are then compared. The chromosome with the highest ranking 

array value is retained, being placed in a temporary holding array, G ( . In case of ties the chromosome 
selected first is retained and placed into the temporary holding array. This process continues until N c 
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chromosomes have been selected. The chromosomes represented by G f are next used by succeeding 
modification operators to produce G n+1 . 

Bin selection — The bin selection algorithm is different from greedy and tournament selection, as 
implemented here, in two general ways. First, the bin selection algorithm chooses its chromosomes from 
the accumulation file, not the active file. Second, the bin selection algorithm divides the distance along 
the pareto front into equal segments or “bins” using an arc-length computation and then selects an equal 
number of chromosomes from each bin using a random process. As with the other options, the selected 
chromosomes are placed in a temporary holding array, G f . 

Since this particular bin selection algorithm requires the computation of arc-length along the pareto front, 
it inherently requires optimization problems with only two objectives, i.e. , the pareto fronts must be curves 
in two dimensions. Another option, which works for a general number of objectives, is available, but has 
not been utilized in the present study. 

The present bin selection algorithm attempts to improve the selection process in two ways. The equal 
partition of the pareto front based on arc-length — not on population — forces increased emphasis for 
selection from underdeveloped regions of the pareto front and reduced emphasis for selection from 
overdeveloped regions. This keeps unwanted clustering from occurring and allows details from all 
portions of the pareto front to develop. Such a selection algorithm might be particularly useful for 
problems that have disjoint or discontinuous pareto fronts, which are often associated with design spaces 
that contain one or more local optima, whether the design space is smooth or not. 

Since the chromosomes that are selected all possess number one rankings, the chance that a superior 
set of chromosomes will be selected using this option is increased, i.e., the chance that the GA will 
converge more quickly is increased. This later characteristic, while likely, is not assured for all 
optimization problems. This is because some pareto fronts are most optimally approached using 
modification operators that require information from the non-optimal design space interior. Knowing which 
type of situation exists a priori is impossible for most realistic engineering design problems. 

Modification Operators 

P Vector — After the new chromosomes have been selected and placed in the temporary holding array, 
G { , they must be modified using one of several modification operators to obtain the ( n + 1) st generation of 

chromosomes, G n+1 . In the present implementation four modification operators are used — passthrough, 
random average crossover, perturbation mutation and mutation. How many chromosomes are modified 
with each operator is controlled by the P vector, which consists of four parameters — p B , p A , p P , p M . 
Each element of the P vector controls one modification operator. The value of each P vector element 
ranges from 0 to 1.0, and, for consistency, the sum of all four elements must always equal one. A P 
vector equal to 0.1, 0.3, 0.3, 0.3, for example, will cause the first 10 percent of the chromosomes to be 
modified using the passthrough operator, the next 30 percent to be modified using random average 
crossover, the next 30 percent to be modified using perturbation mutation and the last 30 percent to be 
modified using mutation. That is, p B =0.1, p A =0.3 , p P = 0.3 , and p M = 0.3 . 

The passthrough operator is always performed first. After passthrough is complete, the implementation 
order of the remaining operators is immaterial. Once all values of G n+1 have been established, the 
algorithm proceeds to fitness evaluation, ranking and then onto succeeding generations until the 
optimization is sufficiently converged. 

Passthrough — The simplest operator used in the present GA is “passthrough.” As the name implies, a 
certain number of chromosomes with the highest individual fitness values are simply “passed through” to 
the next generation from G f to G n+1 without modification. The number of chromosomes that are passed 
through to the next generation is controlled by the first parameter in the P vector, p B . The passthrough 
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operator is always performed on the first pgA/ c chromosomes in G l . Care must be taken when choosing 
p s and N c so that pgA/ c >N 0 . If this is done, the chromosomes with the highest individual fitness values 
will always be passed through to the next generation, thus guaranteeing that none of the individual 
maximum fitness values will ever drop during GA iteration. 


Random average crossover — The next GA modification operator is called random average crossover and 
is implemented by first selecting two random chromosomes and Xj 2 from G ( . Next, the two selected 
chromosomes are combined on a gene-by-gene basis using the following formula: 


n+'l 


-V 

- 9 \ X i j 1 ' 


+ xj J2 ) / = 1,2, --,A/ g 


( 10 ) 


where x" ; +1 corresponds to the /' th gene in the y' th chromosome associated with G n+1 and xjj^ and x/ j2 

correspond to the /' th genes from the randomly chosen chromosomes X^ and Xy 2 . The number of 

chromosomes modified using the random average crossover operator is determined by the parameter 
p A — the second element in the P vector. 


Perturbation mutation — The next GA modification operator is called perturbation mutation and is 
implemented by first selecting a random chromosome X) from G ; . Next, a probability test is performed 

for each gene xjj in the selected chromosome involving a call to a random number generator R(0,1). If 
the returned random number is less than p-i the gene is modified using 

x,r = xjj + (x maX) - x min; )[R( 0, 1)- 0.5 ]p (11) 


where p is a user-specified tolerance that controls the size of the perturbation mutation, and p-i is a user- 
specified control parameter that statistically controls the number of genes that are modified. For sensible 
results the values of p and p-i must be between 0 and 1.0. 


Because this operator can cause the value of a particular gene to exceed one of its constraints (x max . or 
x min . ), checks are required to prevent this. The number of chromosomes modified using the perturbation 
mutation operator is determined by the parameter p P — the third element in the P vector. 


Mutation — The last GA modification operator used in the present study is called mutation and is 
implemented similarly to the perturbation mutation operator. First, a random chromosome Xy is chosen 

from G ( . Next, a probability test is performed for each gene xjj in the selected chromosome involving a 
call to a random number generator R(0,1). If the returned random number is less than p 2 the gene is 
given a completely different value using 

x iJ = (-*max, — )^(0' "I) +-^min, (12) 

The parameter p 2 is a user-specified control parameter that statistically controls the number of genes 
that are modified. For sensible results p 2 must be between 0.0 and 1.0. The number of chromosomes 
modified using the mutation operator is determined by the parameter p M — the fourth element in the P 
vector. 

Gene-space transformation — An option for accelerating GA convergence for multi-objective optimization 
problems is described in this section. When modifying G l to obtain G n+1 it is sometimes advantageous to 
first transform or reorient each chromosome using a gene-space transformation procedure. The 
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modification operators are then applied to the transformed gene values in each chromosome. Once the 
modifications are complete the chromosomes are transformed back using the inverse of the original 
transformation so as to preserve the identity of each gene. 

This gene-space transformation procedure, which can be viewed as a gene-space rotation of 
coordinates, causes a linear coupling between each of the genes, and thus, affects how they are 
changed in the modification process. In some cases the gene-space transformation procedure can 
significantly improve GA convergence. The transformation procedure, which theoretically works for any 
number of objectives, will now be presented for the two-objective case, i.e. , for N 0 = 2. 

The idea behind the transformation procedure is to perform a rotation of coordinates in gene space using 
an angle-preserving, length-preserving orthogonal transformation. For this purpose a simplified version of 
the Gram-Schmidt orthogonalization is used. 42 The current set of r? th generation chromosomes (those 
that have been newly selected and placed in the temporary holding array G f ) can be written as 

*12 *i/v c 1 


x n g ,n c J 

where each column is one of the chromosomes in the active file holding array. It is this matrix that is to be 
transformed using 


*11 

*2,1 

*/V„,1 


G n =UG f (13) 

where U is a unitary matrix of rank N G that needs to be constructed and G n is the resulting transformed 
matrix. To construct U, a set of basis vectors R needs to be established. The elements of R are defined 
by 


// J =X 2 -X 1 for v = 1 


(14a) 


and 


_ J 1 if i = j 
r ‘ J jo if i*j 


for j > 2 


(14b) 


The simple choices made for rjj when j >2 serve to simplify the transformation matrix construction 
without sacrificing overall generality. 


As can be seen from Eq. (14a), the first coordinate direction in the new transformed coordinate system is 
chosen to be parallel to X 2 - X-,. As was mentioned earlier, the chromosome with the best fitness for the 
first objective is always placed into X 1( and the chromosome with the best fitness for the second objective 
is always placed into X 2 . This convention is crucial for the success of the transformation algorithm, as it 
causes the first transformation coordinate to be aligned with the pareto front endpoints. It is imperative 
that the first element of the P vector — the element that controls passthrough — is large enough so that 
both X-i and X 2 are always passed through to the next generation without modification. 
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The unitary transformation matrix U is constructed column by column using 


first column 


second column 


Ui, i = 



Vi = r i,2 -u^ViA . 



n th column 


Vi =r-,n 


n - 1 

jU i. j i 


y=i 



(15a) 


(15b) 


(15c) 


Once G n is obtained the modification operators are applied the same as without transformation to obtain 
G n+1 . The final G n+1 values are then obtained using an inverse transformation indicated by 

G n+1 =U 7 G n+1 (16) 

Because U is a unitary matrix, its inverse is simply its transpose, and thus, it is quite easily constructed. 

The gene space transformation option is controlled using the ITRAN control parameter. If ITRAN=0, no 
gene space transformation in implemented. If ITRAN=1 the gene space transformation algorithm is 
activated. 


Model Problem Descriptions 

To study the relative merits of different GA variations it is useful to use a number of analytic model 
problems. The simplicity of these models — or more appropriately, the speed with which one function 
evaluation can be performed — allows for the removal of statistical variation by averaging a large number 
of tests. The simplicity also allows for the evaluation and comparison of many GA attributes using a 
common basis. In the present study, four multi-objective model problems are used — each with different 
characteristics so as to stress different attributes of any GA algorithm. 

During the discussion of results, design space attributes such as “volume” will be mentioned. For design 
spaces with many dimensions, i.e., many genes, the concept of volume is not a precise one — “hyper- 
volume” being more appropriate. Even the concept of a “hill” in a design space with many dimensions is 
difficult to consider. In the present study terms such as “hill,” “peak” or “volume” will be retained with the 
understanding that the “hyper-” counterparts are, in most cases, more appropriate. 

The first two model problems — MP1 and MP2 — utilize the same formula, introduced in Holst and 
Pulliam. 10 It is based on a simple multi-dimensional paraboloid-like function, which can be replicated with 
different positions and different altitudes to create a large number of multi-modal design spaces. In 
addition, different groupings of this basic function can be utilized to create a large number of multi- 
objective optimization problems. The basic analytic function is given by 
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k=\2,-,N 0 


(17) 


n g 

a m,k =Z (*/ ~ C i,m,kY 

i = 1 

-a„ 


m = 1,2, - • A/, 


b m , k =h, 


1 m.k c 


a m ,k 

Nq 


M 


z k =max(bi' k ,— ,b Nk " k ) J 

where the z ’s are the hill altitudes (design space objectives that need to be maximized), the x ’s are the 
genes (design space decision variables — note that the j subscript has been omitted), the c’s are free 
parameters, the h ’s are peak values for each hill or mode and the a, b quantities are intermediate results. 
The c and h parameter values can either be user input or specified via a random number generator. Once 
established for a particular problem, they do not change. The / subscript is the gene number and varies 
from 1 to N g , the maximum number of genes. The m subscript is the mode number and varies from 1 to 
N m , the maximum number of modes or hills in each design space. Finally, the k subscript is the 
objective number and varies from 1 to A/ 0 , the maximum number of objectives. For the present study, 
N 0 = 2, i.e., only problems with two objectives will be considered. In Eq. (17), the goal is to find values of 
the x ’s that will maximize the z values, and, of course, to do so without using any knowledge one might 
obtain by looking at Eq. (17). With the hill-climbing model presented in Eq. (17), the effect of N G , N M and 
N 0 can be studied, either collectively or individually. 

For simplicity in the present study, the gene limits for each gene for MP1 and MP2 are forced to be the 
same no matter what value of N G is used, i.e., 

■^max 1 = X max 2 
^min, — x min 2 


- x max WQ -2.5 

= ••• =X min„ G =-2.f 


Model Problem No. 1 — The first model problem — MP1 — utilizes Equation (17) with N M =1 and N Q = 2. 
This produces a model hill-climbing problem with a single peak or mode in each objective’s design space. 
The c i m k and h mk coefficients within Eq. (17) are defined as 

C/ii = 0.5, Cj ^ 2 = — 0.5, / = 1,- • A/q 

and 

1 ^= 100 . 0 , = 100.0 

Contour plots showing the shape of MPI’s design space are displayed in Fig. 2. The contour levels of the 
first and second objectives, z-, andz 2 , are plotted in Figs. 2a and 2b, respectively, each as a function of 
the two gene variables, x-i and x 2 . Note that z-, and z 2 are maximized at a position (x-|,x 2 )= (0.5, 0.5) 
and (x 1 ,x 2 )=(-0.5,-0.5), respectively, in keeping with the c imk parameters defined above. 
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Fig. 2 Contours plotted in gene space for a typical two-objective optimization scenario from MP1 , N G = 2 , 
Nm = 1 , A/ 0 = 2 . 


Model Problem No. 2 — The second model problem — MP2 — utilizes Equation (17) with N M =32 and 
N 0 = 2 . This produces a model hill-climbing problem with 32 peaks or modes in each objective’s design 
space. The c j m k and h mk coefficients within Eq. (17) are defined as 


and 


c /,i,i = °-5 j 

c /,i,2 = - 0-5j 
= ^(0.1)(-Xmax, — 0.5) + 0.5 1 
C/, m , 2 = -R(0,1)(x maXj - 0.5) + 0.5 j 


i=\-,N G 

i= \-,N G ,m = 2 ,---,N m 


= 100.0 
\ 2 = 100.0 
Vi = 90.01 
hm,2 = 90.0j 


m = 2, - • • , A/, 


M 


In the above, each time the random number generator R(0,1) is encountered, a new random number is 
generated. With this specification of parameters, the N M -1 local extrema for each objective’s design 
space are randomly placed within the region x min <x ( - <x maXj but excluded from the region 
-0.5 <Xj <0.5. This forces the MP2 pareto front to be the same as the MP1 pareto front. The final 
answer is the same, but convergence to the final answer for MP2 is likely to be more difficult, having to 
traverse a design space with large numbers of secondary peaks. 


Model Problem No. 3 — The third model problem (MP3), introduced in Tanaka, et al., 43 is a minimization 
problem inherently limited to two objectives and two genes. Its analytic specification is given by 


minimize (z-|,z 2 ) 


(18a) 
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where 


*i = *i . z 2 = x 2 


(18b) 


In addition to the above, the following side constraints must be simultaneously satisfied: 


0<x-| ,x 2 <n (19a) 

0 > -x-, 2 -x| + 1 + O.lcos^l 6tan 1 (19b) 

0.5 > (x-|- 0.5f + (x 2 -0.5) 2 (19c) 

Clearly, the complexity of this problem is contained within the constrains, especially the middle constraint 
given by Eq. (19b). 

The contour levels of the two objectives from MP3, z-, andz 2 , are plotted in Figs. 3a and 3b, respectively, 
each as a function of the two gene variables, x-i and x 2 . The lower constraint boundary in each of these 
plots essentially corresponds to the pareto front and is defined by the middle constraint [Eq. (19b)]. The 
upper boundary corresponds to the last constraint [Eq. (19c)]. The linear behavior of each function within 
the constrained feasible design space is obvious. The first objective function z is minimized at the 
extreme left hand side of the feasible domain, and the second objective function z 2 is minimized at the 
extreme bottom. 



a) Objective 1 (z-i) b) Objective 2 (z 2 ) 

Fig. 3 Contours plotted in gene space for the two-objective optimization problem MP3, N G = 2 , N 0 = 2. 


Model Problem No. 4 — The fourth model problem (MP4), introduced in Kursawe, 44 is a minimization 
problem inherently limited to two objectives, N 0 = 2 , but scalable to any number of genes, N G >2. Its 
analytic specification is given by 


minimize (z 1 ,z 2 ) 


(20a) 
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where 



-0.2 

5sinx, 3 ) 


(20b) 


In addition to the above, the following constraints must be simultaneously satisfied: 

-5 < x, < 5 , i = \-,N G (20c) 

Contour plots showing the shape of MP4’s design space for a two-gene, two-objective optimization 
problem are displayed in Fig. 4. The contour levels of the first and second objectives, z-, and z 2 , are 
plotted in Figs. 4a and 4b, respectively, each as a function of the two gene variables, x-, and x 2 . Note 
that the z-, design space is relatively well behaved, similar to an inverted cone in shape. Its minimum 
value (approximately -10.0) occurs at (x 1; x 2 )= (0.0, 0.0) . The z 2 design space is filled with noise and 
represents a nightmare for many optimization methods. Its minimum value (approximately -7.7515) 
occurs at (x.,,x 2 )= (-1.1527,-1.1527) . 

The pareto front for MP4 with two genes is plotted — both in objective space and gene space — in Fig. 5. 
As can be seen, the pareto front for this optimization problem is disjointed and convoluted. The curves 
displayed in Fig. 5 have been determined numerically — from a GA computation with a large number of 
generations — and thus are approximate, but nevertheless, are useful in displaying the complex nature of 
the MP4’s design space. Note how the middle branch of the pareto front plotted in objective space is 
actually split between two locations in gene space. Figure 6 shows additional pareto fronts from MP4 in 
objective space for several values of N G . For each additional gene the pareto front is shifted ten units to 
the left in z 2 and an additional discontinuous branch appears. 




b) Objective 2 (z 2 ) 


Fig. 4 Contours plotted in gene space for a typical two-objective optimization scenario from MP4, N G = 2, 
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a) Objective space b) Gene space 

Fig. 5 Pareto front for MP4 in objective space and gene space, N G =2, N 0 = 2. The red (•), green (□) 
and blue (a) segments in objective space map directly to the red, green and blue segments in gene 
space, respectively. 



Fig. 6 Pareto fronts from MP4 in objective space for several values of N G , N 0 = 2. 


Computed Results 

Area Error Norm for Two-Objective Problems 

In order to evaluate the accuracy and efficiency of an optimization algorithm it is important to have a 
proper error norm to assess the level of convergence. For single objective problems a suitable error norm 
can be obtained using the simple arithmetic difference between the current “best fitness” and the exact 
answer. This, of course, works quite well when working with model problems for which the exact answer 
is known. 

For multi-objective optimization a workable error norm is more elusive to obtain since a range of solution 
values — the so-called pareto front — is being sought. The topic of how to define easy-to-implement and as 
well accurate error norms for comparing one pareto front approximation to another for the same problem 
is discussed at length in Zitzler, et al. and Knowles and Corne. 17 A suitable norm encompassing the 
optimal values of each of the individual objectives is one possibility. However, such a norm would not 
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take into account the trade-off regions of the pareto front. Another possibility is the attainment surface 
approach of Fonseca and Fleming. 45 In this approach a set of equally spaced sampling lines that intersect 
the full breadth of the pareto front are used. Statistical measures of goodness can then be developed 
based upon how many intersection points one pareto front has that are superior to another pareto front. 

In the present study, the area between the current pareto front and the exact pareto front — or an 
accurately computed numerical representation of the exact pareto front — is used as the error norm to 
determine level of convergence. As the size of this quantity goes to zero, the current solution approaches 
the exact solution. Of course, this error norm is only valid for multi-objective problems involving two 
objectives, as a volume computation would be required for three objectives and a hyper-volume 
computation would be required for four or more objectives. As such this report is restricted to multi- 
objective problems in which N 0 = 2 . 

The area is computed numerically by dividing the region between the discretely evaluated exact pareto 
front and the current approximation into a series of triangles as shown in Fig. 7. The area error norm is 
the sum of each of the individual triangular areas. Thus, an approximation to the exact pareto front that 
fails to match the exact solution in any location will produce a non-zero value for the area error norm. As 
the approximate pareto front approaches the exact solution the area error norm approaches zero. 



Fig. 7 Area error norm computational process for a two-objective pareto front. 

As the area error norm approaches zero, the numerical error associated with the triangular discretization 
will eventually become significant. At this point optimization algorithm convergence can no longer be 
monitored using the area error norm. The point at which this occurs is controlled by how many points the 
optimization algorithms finds on the approximate pareto front, which is difficult to control, and by the 
number of points used for the discrete evaluation of the exact pareto front, which is easily controlled. 

To gain an idea about the importance of this truncation error, at least the part associated with exact 
solution discretization, several convergence histories from MP1 using different numbers of points to 
discretize the exact solution (NEXACT) are presented in Fig. 8. Each of these convergence histories is 
averaged over 20 solutions to remove statistical variation (NSEED=20). As can be seen, the effect of 
discretization error on the convergence history curve creeps up to higher and higher values of error as 
NEXACT is reduced. However, even when NEXACT is at a relatively coarse value, e.g., NEXACT=51, 
the present area error norm is able to track the convergence history curve established using larger values 
of NEXACT for nearly four orders of magnitude. Each convergence history curve presented in this report 
will use an NEXACT >1000. 


19 





Fig. 8 Effect of area error norm resolution (NEXACT) on GA convergence, MP1, NSEED = 20, 
ISELECT = 3 , ITRAN = 1, Afe=8, N u = 1, N 0 = 2, N c =20, ,0 = 0.2, p, = 0.5, p 2 = 0.5 and 
P = (0. 1,0.3, 0.3, 0.3). 


Removal of Stochastic Characteristics from Computed Results 

Genetic algorithms are stochastically-based search algorithms and, as such, produce results with 
statistical variation from case to case, even if the only quantity varied is the initializing seed in the random 
number generator. An example of this is displayed in Fig. 9 where three GA convergence histories — area 
error norm versus number of function evaluations — are compared. The number of function evaluations is 
used as a measure of computational work throughout this study despite not including the computational 
work associated with GA algorithm overhead, because it is easy to define and because it does not vary 
from computer to computer. For most applications, the computational work associated with the GA 
optimization overhead is easily dominated by the function evaluation aspect of the computation, and thus, 
the present results are useful in determining which GA parameter and design space attributes produce 
the most efficient computational results. 

All three of the convergence histories presented in Fig. 9 are from MP3. The GA parameter values are 
given by ISELECT = 1, ITRAN = 0 , N G =3, N 0 = 2, A/ c =20, ,0 = 0.2, p,=0.5, p 2 = 0.5 and 
P = (0.1, 0.3,0.3, 0.3). This set of GA parameter values defines the so-called baseline solution for this 
study. The GA parameters utilized in this solution — all except ISELECT and ITRAN — were determined via 
trial and error and represent the most efficient set of GA parameters for this model problem. All GA 
parameters not being varied in this report utilize the default values that are established above. 

Except for the seed value used to initialize the random number generator at the beginning of each GA 
solution, all GA and problem parameters are the same for the three convergence histories displayed in 
Fig. 9. As can be seen, the three convergence histories are different — in some locations the differences 
are significant. These differences are caused by statistical variations encountered during solution 
execution and are typical for a GA search process. For multi-modal cases or cases with noise in the 
design space, the differences can be much larger than those displayed in Fig. 9. 
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Fig. 9 Three GA convergence histories utilizing different initial seed values with all other GA and design 
space parameters held fixed, MP3, ISELECT = 1 , ITRAN = 0 , N G = 3 , N M = 1, N 0 = 2, N c =20, /? = 0.2 , 
Pl =0.5, p 2 = 0.5 and P =(0.1, 0.3, 0.3, 0.3). 

To study the relative effects of various GA or design space parameters on GA convergence, it is 
important to remove this statistical variation — at least most of it — so that the average effect of each 
parameter being studied can be accurately ascertained. This is accomplished by running each case 
many times with different initializing seed values, and then averaging the results. The number of 
independent solutions that must be included in the average to produce results free of statistical 
variation — NSEED — can be determined by computational experiment. Results for such an experiment are 
presented in Fig. 10. 



Fig. 10 Pareto front convergence histories averaged over a different number of solutions (NSEED 
values), all other GA and design space parameters held fixed, MP3, ISELECT = 1, ITRAN = 0, N G =3, 
N m - 1 , N 0 =2, N c =20, /? = 0.2, Pl =0.5, p 2 = 0.5 and P =(040.3,0.3,0.3). 

As can be seen, sample-size independence is obtained for an NSEED value in the range of 10-20 for this 
set of conditions. As such, all results presenting convergence efficiency information will be averaged over 
20 solutions, i.e., NSEED = 20 . It should be pointed out that this result is only valid for smooth or single- 
mode design spaces, i.e., design spaces that have only a single optimum and that are not convoluted. 
Noisy or multi-mode design spaces require even larger NSEED values for sample-size independence. 
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Convergence Efficiency Comparisons 


The effect of the selection algorithm (ISELECT) and the gene space transformation procedure (ITRAN) 
on GA convergence efficiency is presented in Figs. 11-13 for the first model problem (MP1). Results for 
three different values of ISELECT are presented: ISELECT=1, greedy selection, ISELECT=2, tournament 
selection, and ISELECT=3, bin selection. In addition, results for the untransformed GA scheme 
(ITRAN=0) and the transformed GA scheme (ITRAN=1) are presented. All convergence history curves 
are averaged over 20 cases, i.e. , NSEED = 20. 

Results showing the effect of two different error norms on the GA convergence performance are 
displayed in Fig. 11. Figure 11a shows convergence histories utilizing the above described area error 
norm plotted versus the number of function evaluations. Figurellb shows convergence histories utilizing 
the endpoint error norm, which is derived solely from the pareto front endpoints and is defined as 

endpoint_error = Z| z M^act ‘ Z ™ X \ 

k= 1 


where z™* act and z™ ax are the pareto front endpoints associated with the exact and approximate 
solutions, respectively. Keep in mind that the solutions used to form the average convergence history 
curves displayed in Figs. 11a and 11b are the same. The differences that are seen in these two different 
sets of curves are solely associated with how the error norms were computed. 

As can be seen from Fig. 11, conclusions about which scheme variations are the most efficient are 
different based upon which error norm is used. Because the endpoint error norm does not take into 
account the entire pareto front, it is grossly inaccurate for many scheme variations, showing better than 
expected convergence efficiency for selection schemes which favor pareto endpoints and poorer than 
expected convergence efficiency for selection schemes that do not favor pareto endpoints. 



a) Area error norm 



b) Endpoint error norm 


Fig. 11 Pareto front convergence histories averaged over different ISLECT and ITRAN values, all other 
GA and design space parameters held fixed, MP1, NSEED = 20, N G =8, N M ='\, N a = 2, N c =20, 
>3 = 0.1, pi =0.5, p 2 = 0.1 and P = (0. 1,0.3, 0.3, 0.3). 
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a) 10 generations 


a) 10 generations 




b) 100 generations b) 100 generations 



c) 1000 generations 


Fig. 12 Solution development as seen from 
objective space for MP1, ISELECT = 1, 
ITRAN = 0 , Ne=8, N m = 1, N 0 = 2, N c =20, 
P = 0.2, p, =0.5, p 2 = 0.5 , P =(0.1, 0.3, 0.3, 0.3). 



c) 1000 generations 

Fig. 13 Solution development as seen from 
objective space for MP1, ISELECT = 3, 
ITRAN = 1, N G = 8, N m = 1, N 0 =2, N c =20, 
p= 0.2, p, =0.5, p 2 = 0.5 , P =(0.1, 0.3, 0.3, 0.3). 
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Confirmation of this observation can be obtained by looking at the development of the pareto front at 
different GA generation levels. This is done for two cases in Figs. 12 and 13. Figure 12 shows pareto 
front development at the 10, 100 and 1000 generation levels for the ISELECT = 1, ITRAN = 0 case. 
Figure 1 3 shows pareto front development at the same generation levels for the ISELECT = 3 , ITRAN = 1 
case. Three sets of data are displayed in each plot: the exact pareto front as a solid red line, the 
approximate pareto front as a series of blue circles and the design space — all solutions explored to that 
point — as a series of small green dots. 

As can be seen from Figs. 12 and 13, the ISELECT = 3, ITRAN = 1 case develops its pareto front much 
more rapidly than the ISELECT = 1, ITRAN = 0 case. In fact, the approximate pareto front at 100 
generations for the ISELECT = 3 , ITRAN = 1 case (Fig. 13b) more accurately represents the exact pareto 
front than the approximate pareto front at 1000 generations for the ISELECT = 1, ITRAN = 0 case (Fig. 
12c). This behavior is in general agreement with the convergence history results that are displayed in Fig. 
11, which utilized the area error norm. 

From Fig. 11 it can be seen that the three convergence history curves that utilize the transformation 
procedure (ITRAN = 1) display better convergence performance than the three convergence history 
curves that do not use the transformation procedure (ITRAN = 0). In particular, after an error drop of two 
orders of magnitude, the difference in the number of function evaluations is more than an order of 
magnitude. More importantly, the slope on each of the ITRAN = 1 convergence history curves is steeper 
than any of the ITRAN = 0 curves, indicating that the method speed up will only increase as the GA 
process converges further. 

The effect of the selection algorithm (ISELECT) and the gene space transformation procedure (ITRAN) 
on GA convergence efficiency for MP2 is presented in Fig. 14. Results for all three selection algorithms 
and both values of ITRAN are presented. In each case the pareto front error using the area error norm is 
plotted against the number of function evaluations. Because of the effective increase in design space 
noise caused by using N M =32, all convergence history curves were averaged over 100 solutions, i.e., 
NSEED = 100 . 



Fig. 14 Pareto front convergence histories for 
different ISLECT and ITRAN values, all other 
GA and design space parameters held fixed, 
MP2, NSEED = 100 , N G =8, N M = 32 , N 0 = 2, 
N c =20, >0 = 0.1, pi =0.5, p 2 = 0.1 and 

P = (0.1, 0.3, 0.3, 0.3). 



Fig. 15 Pareto front convergence histories for 
different ISLECT, ITRAN and N M values, all 
other GA and design space parameters held 
fixed, MP1 and MP2, NSEED = 100, Afe=8, 
N 0 = 2 , N c =20, /? = 0.1, Pi = 0.5 , p 2 = 0.1 and 
P = (0.1, 0.3, 0.3, 0.3). 
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For this case the results that utilize bin selection are clearly superior to the results that utilize either 
greedy selection or tournament selection. In addition, the ITRAN = 1 results when used in conjunction with 
bin selection are more efficient. 

Selected convergence history results from Fig. 11a (rerun with NSEED = 100) and Fig. 14 — the fastest 
and slowest curves from each figure — are cross-plotted versus each other in Fig. 15. This allows a 
quantitative look at what the N M parameter does to GA convergence. As can be seen, the faster 
results — ISELECT = 3, ITRAN = 1 — are only slightly inconvenienced by N M =32, whereas the slower 
results — ISELECT = 1, ITRAN = 0 — suffer from the existence of so many local extrema. 

The effect of the selection algorithm (ISELECT) and the gene space transformation procedure (ITRAN) 
on GA convergence efficiency for MP3 is presented in Figs. 16-18. Results for all three selection 
algorithms and both transformation variations are included. For each case presented in Fig. 16, the 
pareto front error using the area error norm is plotted against the number of function evaluations. All 
convergence history curves are averaged over 100 cases, i.e., NSEED = 100. 

As seen from Fig. 16, the ISELECT=3, ITRAN = 1 and ISELECT=3, ITRAN = 0 schemes are still the 
most efficient, but now, there is less variation in convergence efficiency over all of the GA variations 
plotted. Overall the total error drop over 100,000 function evaluations ranged from about one-and-a-half 
orders of magnitude to a little less than two orders of magnitude. For MP1 and MP2 the total error drop 
over the same number of function evaluations ranged from about one to about four orders of magnitude. 
Generally, the convergence efficiency performance for the slower schemes across MP1, MP2 and MP3 
was about the same. The faster schemes for MP1 and MP2, however, suffered dramatically when they 
were applied to MP3. The reason for this is most likely associated with MP3’s nonlinear side constraint, 
which produces a discontinuous and mildly convoluted pareto front. 



NUMBER OF FUNCTION EVALUATIONS 


Fig. 16 Pareto front convergence histories for different ISLECT and ITRAN values, all other GA and 
design space parameters held fixed, MP3, NSEED = 100, N G =2, A/ 0 = 2 , N c =20, /3 = 0.2, p^= 0.5, 
p 2 = 0.1 and P=(0. 1,0.3, 0.3, 0.3). 


Comparisons between the computed and the exact pareto fronts superimposed on the numerically 
generated design space for MP3 are presented in Figs. 17 and 18. Two different cases, corresponding to 
ISELECT = 3 , ITRAN = 0 and ISELECT = 3 , ITRAN = 1, are included. In both cases the computed pareto 
front is represented with red circular symbols and the exact pareto front is represented with green dots. 
The design space data — which includes all solutions that were evaluated during the course of GA 
iteration and that did not violate one of the side constraints — are plotted using a series of black dots. 
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Note that both cases produce good approximations to the exact pareto front. The design space for the 
ITRAN = 0 case is approximately laid out within a triangle — one side along the pareto front, the other two 
along the x-, = constant and x 2 = constant directions. The ITRAN = 1 case design space is more so 
flattened against the pareto front. This is a direct result of the gene-space transformation, as the 45° line 
that stretches from the pareto front endpoints becomes the preferred direction. 
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0 0.2 0.4 0.6 0.8 1 1.2 

OBJECTIVE 1 

b) Computed pareto front in relation to the 
computed design space 

Fig. 17 Pareto front comparisons for MP3, 2500 
generations, ISELECT = 3, ITRAN = 0, A/ G = 2 , 
N 0 = 2, N c =20, p = 0.1, pi =0.5, p 2 = 0.1 and 
P =(0.1, 0.3, 0.3, 0.3). 





0 0.2 0.4 0.6 0.8 1 1.2 

OBJECTIVE 1 

b) Computed pareto front in relation to the 
computed design space 

Fig. 18 Pareto front comparisons for MP3, 2500 
generations, ISELECT = 3, ITRAN = 1, A/ G =2, 
N 0 = 2, N c =20, p = 0.1 , pi =0.5, p 2 = 0.1 and 
P =(0.1, 0.3, 0.3, 0.3). 
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The effect of the selection algorithm (ISELECT) and the gene space transformation procedure (ITRAN) 
on GA convergence efficiency for MP4 is presented in Figs. 19-20. For each case presented in Fig. 19, 
the pareto front error using the area error norm is plotted against the number of function evaluations. All 
convergence history curves are averaged over 100 cases, i.e., NSEED = 100. Results for all three 
selection algorithms and both transformation variations are included for the two-gene variation of this 
model problem. Because of the complication of this model problem, the so-called “exact” solution used in 
the area error norm computation is actually a tightly converged numerical solution run for 10 7 function 
evaluations. 

As seen from Fig. 19, the three convergence history curves that correspond to the ITRAN = 1 results are 
the slowest and the three untransformed curves are the fastest. Because of the convoluted nature of this 
model problem’s pareto front, the transformation procedure actually harms GA convergence. For this 
case the ISELECT = 3, ITRAN = 0 case is the most efficient, achieving a three-and-a-half order of 
magnitude reduction in the area error norm over 100,000 function evaluations. 

The solution evolution for this most efficient case (ISELECT=3, ITRAN = 0) is displayed in Fig. 20. 
Separate solutions, plotted in both objective space and in gene space, are presented at 10, 100 and 1000 
generations. Three data sets are displayed at each generation level, the computed pareto front (red 
circular symbols), the “exact” pareto front (green dots) and the design space (black dots). Note that the 
scale for the gene space plots is greatly expanded. Thus, the gene space plots do not display as many 
points as the objective space plots. As can be seen from Fig. 20, development of the pareto front for this 
model problem is more difficult because of its convoluted nature. 



Fig. 19 Pareto front convergence histories for different ISLECT and ITRAN values, all other GA and 
design space parameters held fixed, MP4, NSEED = 100, N G =2, A/ 0 = 2 , N c =20, p = 0.2, =0.5, 

p 2 = 0.1 and P = (0.1, 0.3, 0.3, 0.3). 
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a) Objective space, 10 generations. 



b) Gene space, 10 generations. 




c) Objective space, 100 generations. d) Gene space, 100 generations. 



e) Objective space, 1000 generations. f) Gene space, 1000 generations. 

Fig. 20 Computed pareto fronts in objective and gene space at selected levels of convergence, MP4, 
ISELECT = 3 , ITRAN = 0 , N G =2, N 0 = 2, N c =20, j3= 0.2, p, =0.5, p 2 = 0.1 and P = (0.1, 0.3, 0.3, 0.3) . 
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Effect of modification operator parameter variation on GA convergence — The effect of the various 
modification operator parameters on GA convergence efficiency is studied in this section. In particular, 
three parameters will be considered, /?, p i and p 2 . The /? parameter is used to control the size of 
perturbations in the perturbation mutation operator. The p-i parameter is used to control the probability 
that any given gene will be perturbed in the same operator. The p 2 parameter is used to control the 
probability that any given gene will be mutated using the global mutation modification operator. The 
number of chromosomes modified by these two operators is controlled by the third and fourth elements of 
the P vector, respectively. For all results in this section, P =(0.1, 0.3,03, 0.3) . Thus, the perturbation 
mutation and global mutation operators are each used to update 30% of all chromosomes. 

Figure 21 shows the effect of j3 on GA convergence efficiency for each of the four model problems. In 
each case the pareto front error using the area error norm is plotted against the number of function 
evaluations. For all results in this section the bin selection algorithm without gene-space transformation is 
used, i.e. , ISELECT = 3 , ITRAN = 0 . All convergence history curves are averaged over 100 cases. 




a) MP1, N g =8. 


c) MP3, Nq =2. 




b) MP2, N G = 8 . d) MP4, N G = 3 . 


Fig. 21 Pareto front convergence histories for different model problems showing the effect of /? variation 
on GA convergence efficiency, ISELECT = 3, ITRAN = 0, NSEED = 100, N 0 = 2, N c =20, p 1 = 0.5, 
p 2 = 0.1 and P =(0.1, 0.3, 0.3, 0.3). 
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As can be seen from in Fig. 21 , the convergence history curves associated with different values of p tend 
to cross each other. This trend exists for each model problem and is manifested by the fact that the 
optimal value of p varies as a function of convergence level. But the trend for which values are most 
efficient and which are least efficient changes across the various model problems. For MP1 and MP2, 
which have a simple, well-behaved pareto front, moderate values of p are most efficient for early 
convergence. As the final pareto front is approached, smaller and smaller perturbations are required for 
optimal convergence. This is the same trend observed for single-objective optimization problems as 
reported in Ref. 10. 

For MP3 and MP4, which have more complex pareto fronts, being either mildly convoluted or strongly 
convoluted, the smaller values of p are generally best for early convergence, and the larger values of p 
are generally best for late convergence. The precise reason for this trend is unclear, but it suggests that a 
GA approach that uses adjustable values of p might be attractive. 

Figure 22 shows the effect of p- ( on GA convergence efficiency for each of the four model problems. In 
each case the pareto front error using the area error norm is plotted against the number of function 
evaluations. For all results in this section the bin selection algorithm without gene-space transformation is 
used, i.e. , ISELECT = 3 , ITRAN = 0 , and all convergence history curves are averaged over 100 cases. 




b) MP1, N g =8. 


c) MP3, Ng =2 




b) MP2, Nq = 8. 


d) MP4, Nq =3 


Fig. 22 Pareto front convergence histories for different model problems showing the effect of p-i variation 
on GA convergence efficiency, ISELECT = 3, ITRAN = 0, NSEED = 100, A/ 0 = 2 , A/ c = 20 , p= 0.2, 
p 2 = 0.1 and P = (0.1, 0.3, 0.3, 0.3). 
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As can be seen from Fig. 22, the effect of p-, on GA convergence efficiency is small for MP1 , MP2 and 
MP4 and only marginal for MP3. For MP1, MP2 and MP4, intermediate values of p-i are generally best, 
while for MP3 larger values are best. 

Figure 23 shows the effect of p 2 on GA convergence efficiency on two of the model problems — MP1 and 
MP4 — the two model problems with the simplest pareto front and the most complex pareto front. In each 
case the pareto front error using the area error norm is plotted against the number of function 
evaluations. For each plot, results for ITRAN = 0 and ITRAN = 1 are presented for just two values of p 2 — 
0.0 and 1.0. 

As can be seen in Fig. 23, the effect of p 2 on GA convergence efficiency is negligible for all cases 
presented. Many other cases were run, which utilized different values of p 2 . The effect of p 2 on GA 
convergence efficiency was always negligible. This is in direct contrast with the single-objective results 
obtained in Ref. 10. For the single-objective results, the effect of p 2 on GA convergence efficiency was 
also negligible for cases in which design space noise was nonexistent or small, but dramatic for cases in 
which the design space noise was significant. It can only be surmised that none of the cases attempted in 
this study possessed enough design space noise to demonstrate the effect of p 2 on GA convergence 
efficiency. 
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a) MP1, N g =8. 


b) MP4, N g =3. 


Fig. 23 Pareto front convergence histories for different model problems showing the effect of p 2 variation 
on GA convergence efficiency, ISELECT = 3, NSEED = 100, N 0 = 2, N c =20, 0 = 0.2, pi=0.5 and 
P = (0.1, 0.3, 0.3, 0.3). 


Effect of P vector variation on GA convergence — The effect of the P vector, which controls which 
modification operators are to be used in going from G n to G n+1 , is studied in this section. Selected 
convergence history results for each of the model problems (MP1; MP2, MP3 and MP4) are displayed in 
Fig. 24. In each case convergence history results for the baseline value of the P vector, 
P=(0. 1,0.3, 0.3, 0.3), along with an assortment of other P vector values designed to systematically 
emphasize each of the later three modification operators — random average crossover, perturbation 
mutation and global mutation — are displayed. For all results in this section the bin selection algorithm 
without gene-space transformation is used, i.e., ISELECT = 3 , ITRAN = 0 . All convergence history curves 
are averaged over 100 cases, i.e., NSEED = 100 . 
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As can be seen from Fig. 24 the effect of different P vector values on GA convergence is not as large as 
for some of the other parameters studied. The variation in efficiency between the best and worst 
convergence history curves — as measured by the ratio of the number of function evaluations at fixed 
error — is generally between three and five for each of the model problems. In addition, the asymptotic 
convergence rates for each P vector for each model problem are almost identical, thus indicating that the 
relative efficiency comparisons will continue to smaller values of error. 

Additional trends across each of the model problems are also observable. The two best convergence 
history curves for each model problem are always associated with the baseline P vector and the P 
vector that emphasizes the perturbation mutation operator — P = (0.1, 0.1, 0.7, 0.1) . The P vector that 
emphasizes global mutation — P =(0.1, 0.1, 0.1, 0.7) — is always associated with the worst convergence 
history efficiency. 




a) MP1, Nq =8. 


c) MP3. 




b) MP2, Nq =8. 


d) MP4, N g =3. 


Fig. 24 Pareto front convergence histories for different model problems showing the effect of P vector 
variation on GA convergence efficiency, ISELECT = 3, ITRAN = 0, NSEED=100, N 0 = 2, N c =20, 


J3= 0.2, pi = 0.5 and p 2 = 0.1. 
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Effect of number of genes on GA convergence — The effect of the number of genes used to define the 
design space ( N G ) is studied in this section. Selected convergence history results for two of the model 
problems (MP1 and MP4) are displayed in Fig. 25. For MP1 N G is varied from 2 to 64 and for MP4 N G is 
varied from 2 to 8. For all results in this section the bin selection algorithm without gene-space 
transformation is used, i.e., ISELECT = 3, ITRAN = 0. All convergence history curves are averaged over 
100 cases, i.e., NSEED = 100. 

As can be seen from Fig. 25, the effect of N G on GA convergence is dramatic. The key thing to look at in 
Fig. 25 is the slope of each curve, i.e., how much does each curve drop during the course of 100,000 
function evaluations? Not surprisingly, as N G increases, the rate at which the GA converges decreases. 
Thus, in choosing a design space discretization, it is imperative to select a parameterization that leads to 
the smallest number of decision variables as possible. 
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a) MP1. 


b) MP4. 


Fig. 25 Pareto front convergence histories for different model problems showing the effect of N G variation 
on GA convergence efficiency, ISELECT = 3, ITRAN = 0, NSEED = 100, A/ 0 = 2 , A/ c =20, /?= 0.2, 
Pl =0.5, p 2 = 0.1 and P =(0.1, 0.3, 0.3, 0.3) . 


Effect of number of chromosomes on GA convergence efficiency — For all GA algorithm variations studied 
in this report, the number of chromosomes in each population, N c , is held fixed over the entire 
optimization process. Nevertheless, the value of N c is a user specified parameter that can have an 
impact on GA convergence. Quantification of this impact for the four model problems is presented in Fig. 
26. In each case the error in the pareto front using the area error norm is plotted against the number of 
function evaluations. For all results in this section the bin selection algorithm without gene-space 
transformation is used, i.e., ISELECT = 3, ITRAN = 0. All convergence history curves are averaged over 
100 cases, i.e., NSEED = 100. 

As can be seen from Fig. 26, especially Figs. 26a and b, which are associated with the simple pareto 
fronts of MP1 and MP2, the initialization penalty increases with increasing values of N c . Nevertheless, the 
asymptotic convergence is not affected by the value of N c . In particular, after enough function 
evaluations, somewhere around 5000-10000, all of the convergence history results for a given model 
problem collapse to a single curve. Generally, the smaller values of N c produced the fastest 
convergence. This has a significant implication for massively parallel applications and suggests that GA 
optimization approaches can be parallelized on virtually any number of processors with an extremely high 
degree of efficiency. 
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a) MP1, N g = 8. 


c) MP3, N g = 2. 




b) MP2, N G = 8 . 


d) MP4, N g =3. 


Fig. 26 Pareto front convergence histories for different model problems showing the effect of number of 
chromosomes, N c , on GA convergence efficiency, ISELECT = 3, ITRAN = 0, NSEED = 100, A/ 0 = 2, 
Pl =0.5, p 2 = 0.1 and P =(0.1, 0.3, 0.3, 0.3) . 


Conclusions 


A genetic algorithm (GA) optimization procedure, with several selection algorithm options and a new 
gene-space transformation procedure, has been presented. The new algorithm is especially designed for 
multi-objective optimization problems. It uses real-number encoding to represent all design space 
decision variables as genes and populations of fixed size to go from generation to generation. Four 
modification operators are utilized to advance from one generation to the next. They include passthrough, 
random average crossover, perturbation mutation and mutation. The standard output for this approach is 
a pareto front, which includes the best solutions from each objective, as well as a range of non-dominated 
“tradeoff” solutions in between. 
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The GA optimization procedure converged to the global pareto front optimum for every case attempted, 
demonstrating robustness and wide applicability. In some cases convergence was achieved quickly, 
while in other cases convergence was much slower. Four model problems with different design space 
characteristics were used to systematically study the effects of the various GA parameters on 
convergence performance. The following specific conclusions were obtained: 

Of the three selection algorithms tested, greedy, tournament and binning, the latter produced the most 
efficient convergence across all of the problems studied. The gene-space transformation procedure was 
successful for model problems with simple pareto fronts, achieving improvements in convergence 
efficiency. It was moderately successful for model problems with mildly convoluted pareto fronts and 
displayed a degradation in convergence performance for pareto fronts that were strongly convoluted. 

Key design space features that tended to diminish GA convergence performance were large numbers of 
genes and design spaces that were highly convoluted. The fact that these areas would slow convergence 
is intuitive. A more difficult less intuitive area to assess is the effect of GA parameters on convergence 
performance. The conclusions in this area are presented by looking at each GA parameter separately. 

In general, GA parameter values had only a small to moderate effect on GA convergence efficiency. The 
perturbation mutation parameter controlling size of the perturbation mutations (p) displayed a crossover 
behavior for each model problem studied, with one value of p more efficient during early convergence 
and another value more efficient during late convergence. The perturbation mutation probability (p-,) 
generally had a small effect on GA convergence efficiency with moderate values producing the best 
convergence. The global mutation probability (p 2 ) had a negligible effect on GA convergence efficiency 
for all cases tested. Lastly, the P vector had a moderate effect on convergence efficiency with the default 
value — P = (0.1, 0.3, 0.3, 0.3) — or values that emphasized the perturbation mutation operator being the 
best. 

The number of chromosomes used in each generation, N c , had a small effect on GA convergence 
efficiency for all situations studied. Initially, the larger values of N c seemed to possess an initialization 
penalty that suggested their convergence efficiency would be degraded, but asymptotically, the larger 
values of N c produced the same convergence rate as smaller values of N c for all model problems 
studied. This suggests that the GA optimization procedure would be extremely attractive for 
implementation on massively parallel computers. 
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